Data was simulated using VirtualCommunity code.
Simulated data contains 20 data sets.
NB - if chain/sample parameters is modified - should be changed in the whole section(dataset)
data<-sim_data$EnvEvenSp5
data <- list(
Y = subset(data, select = -env),
X = cbind(1, scale(poly(data$env, 2))),
covx = cov(cbind(1, scale(poly(data$env, 2)))),
K = 3,
J = ncol(data) - 1,
n = nrow(data),
I = diag(ncol(data) - 1),
df = ncol(data)
)
Y_cor<-cor(data$Y)
to_prec<-function(m){
n<-dim(m)[1]
Tau_n<-matrix(nrow=n, ncol=n)
for (j in 1:n) {
for (k in 1:n){
Tau_n[j, k] <- -m[j, k]/sqrt((m[j,j]*m[k,k]))
}
}
return(Tau_n)
}
#Tau_n<-matrix(nrow=dim(model$mean$Tau)[1], ncol=dim(model$mean$Tau)[1])
Tau_n<-to_prec(me5$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))
par(mfrow=c(2,4),oma = c(3, 1, 2, 1))
cols = colorRampPalette(c("dark blue","white","red"))
col2 <- colorRampPalette(c("#4393C3", "#2166AC", "#053061",
"#FDDBC7", "#FFFFFF", "#D1E5F0", "#92C5DE",
"#67001F", "#B2182B", "#D6604D", "#F4A582"))
corrplot(Y_cor, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Correlation cor(Y)")
corrplot(me5$mean$EnvRho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("EnvRho")
corrplot(me5$mean$EnvRho*(!me5$overlap0$EnvRho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("EnvRho signif")
corrplot(me5$mean$Rho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Rho")
corrplot(me5$mean$Rho*(!me5$overlap0$Rho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Rho signif")
corrplot(Tau_n, diag = FALSE, order ="original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Tau")
corrplot(Tau_n*(!me5$overlap0$Tau), diag = FALSE, order ="original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Tau signif")
#j_metric_e5<-metrics_jsdm(me5)
#cat("Success rate ")
#j_metric_e5$success_env
cat(sprintf("Success rate: %s\n", metrics_jsdm(me5)$success_env))
## Success rate: 0.9
##to setup chains parameters
data<-sim_data$EnvEvenSp5
#Rho_sign<-fit_gjam(data,2000,1000,"./gjam_models/gjam5env.rda")
Rho_sign<-load_gjam(data,name="./gjam_models/gjam5env.rda")
#g_metric_e5<-metrics_gjam(Rho_sign)
#cat("Success rate ")
#g_metric_e5$success_env
cat(" ")
cat(sprintf("Success rate: %s\n", metrics_gjam(Rho_sign)$success_env))
## Success rate: 0.8
#setwd("~/Tesi/Code/Ecology-models-master/simcoms-master/ExampleFiles")
fit_hmsc<-function(data,label="Fit",nsamples = 1000,nchains=2,name="./HMmodels/hmtemp.rda" ){
if (label=="Fit"){
Y_data = subset(data, select = -env)
ns<- ncol(Y_data)
np <- nrow(Y_data)
X<-scale(poly(data$env[1:np], 2))
colnames(X)<-c("env","env2")
studyDesign = data.frame(sample = as.factor(1:np))
rL = HmscRandomLevel(units = studyDesign$sample)
m = Hmsc(Y=as.matrix(Y_data), XData=as.data.frame(X), XFormula=~env+env2, distr="probit",
studyDesign = studyDesign, ranLevels = list(sample = rL))
m = sampleMcmc(m, nsamples, thin=10, adaptNf=c(200,200), transient=500,nChains=nchains ,verbose=F)
save(m, file = name)
return(m)
}
if (label=="Load"){
return(load_object(name))
}
}
data<-sim_data$EnvEvenSp5
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm5env.rda" )
#hm_mod<-load_object("./HMmodels/hm5env.rda")
Convergence:
hm_conv<-function(mod){
codaList = convertToCodaObject(mod)
#convergence histograms
hist(effectiveSize(codaList$Beta), main="ess(beta)",lwd=2,col=gray(.6))
hist(gelman.diag(codaList$Beta,multivariate=FALSE)$psrf,lwd=2,col=gray(.6), main="psrf(beta)")
hist(effectiveSize(codaList$Omega[[1]]), main="ess(omega)",lwd=2,col=gray(.6))
hist(gelman.diag(codaList$Omega[[1]], multivariate=FALSE)$psrf,lwd=2,col=gray(.6), main="psrf(omega)")
}
hm_conv(hm_mod)
Study of interactions
metrics_hmsc<-function(Rho, comp=NULL,fac=NULL,only_env=T){
#{
if(!is.null(comp) & !is.null(fac)){
fac_comp <- fac-comp
TP_comp <- FN_comp <-TP_fac <- FN_fac <- FP <- TN <- wrong <- 0
for(i in 1:(nrow(Rho)-1)){
for(j in (i+1):ncol(Rho)){
if(fac_comp[i,j]>0 & Rho[i,j] >0)
TP_fac=TP_fac+1
if(fac_comp[i,j]>0 & Rho[i,j] == 0)
FN_fac=FN_fac+1
if(fac_comp[i,j]==0 & Rho[i,j] == 0)
TN=TN+1
if(fac_comp[i,j]==0 & Rho[i,j] != 0)
FP=FP+1
if(fac_comp[i,j]< 0 & Rho[i,j] < 0)
TP_comp=TP_comp+1
if(fac_comp[i,j]<0 & Rho[i,j] == 0)
FN_comp=FN_comp+1
if((fac_comp[i,j]<0 & Rho[i,j] > 0) | (fac_comp[i,j]>0 & Rho[i,j] < 0))
wrong=wrong+1
}
}
}else{
if(!is.null(comp)){
TP_comp <- FN_comp <- FP <- TN <- wrong <- 0
TP_fac <- FN_fac<-NULL
for(i in 1:(nrow(Rho)-1)){
for(j in (i+1):ncol(Rho)){
if(comp[i,j]>0 & Rho[i,j] <0)
TP_comp=TP_comp+1
if(comp[i,j]>0 & Rho[i,j] >0)
wrong=wrong+1
if(comp[i,j]>0 & Rho[i,j] == 0)
FN_comp=FN_comp+1
if(comp[i,j]==0 & Rho[i,j] == 0)
TN=TN+1
if(comp[i,j]==0 & Rho[i,j] != 0)
FP=FP+1
}
}
# }
# }
}
if(!is.null(fac)){
TP_fac <- FN_fac <- FP <- TN <- wrong <- 0
TP_comp <- FN_comp<-NULL
for(i in 1:(nrow(Rho)-1)){
for(j in (i+1):ncol(Rho)){
if(fac[i,j]>0 & Rho[i,j] >0)
TP_fac=TP_fac+1
if(fac[i,j]>0 & Rho[i,j] <0)
wrong = wrong+1
if(fac[i,j]>0 & Rho[i,j] == 0)
FN_fac=FN_fac+1
if(fac[i,j]==0 & Rho[i,j] == 0)
TN=TN+1
if(fac[i,j]==0 & Rho[i,j] != 0)
FP=FP+1
}
}
}
}
if(only_env){
FP <- TN <- 0
TP_fac <- FN_fac<-TP_comp <- FN_comp<-wrong<-NULL
for(i in 1:(nrow(Rho)-1)){
for(j in (i+1):nrow(Rho)){
if(Rho[i,j] != 0)
FP=FP+1
if(Rho[i,j] == 0)
TN=TN+1
}
}
}
success_env<-TN/(TN+FP)
if(!is.null(TP_comp)){success_comp <- TP_comp/(TP_comp+FN_comp)}else{ success_comp = NULL}
if(!is.null(TP_fac)) {success_fac <- TP_fac/(TP_fac+FN_fac)}else{ success_fac = NULL}
list("FP"=FP,"TN"=TN,"TP_comp"=TP_comp,"FN_comp"=FN_comp,"TP_fac"=TP_fac,"FN_fac"=FN_fac,
"wrong"=wrong,"success_env"=success_env,"success_comp"=success_comp,"success_fac"=success_fac)
}
hm_inter<-function(mod, nchains=2,nsamples = 1000, interact=diag(ns)){
getOmega = function(a,r=1)
return(crossprod(a$Lambda[[r]]))
ns<-mod$ns
postOmega1 = array(unlist(lapply(mod$postList[[1]],getOmega)),c(ns,ns,mod$samples))
postOmega2 = array(unlist(lapply(mod$postList[[2]],getOmega)),c(ns,ns,mod$samples))
postOmega<-abind(postOmega1,postOmega2,along=3)
postOmegaMean = apply(postOmega,c(1,2),mean)
postOmegaUp=apply(postOmega,c(1,2),quantile,0.95)
postOmegaLo=apply(postOmega,c(1,2),quantile,0.05)
postR<-array(dim=c(ns,ns,nchains*nsamples))
for(i in 1:dim(postOmega)[3])
postR[,,i]<-stats::cov2cor(postOmega[,,i])
postRMean = apply(postR,c(1,2),mean)
postRUp=apply(postR,c(1,2),quantile,0.95)
postRLo=apply(postR,c(1,2),quantile,0.05)
Tau = solve(postOmegaMean)
Tau_n = -cov2cor(Tau)
Toplot_R<-postRMean*(!(postRUp>0 & postRLo<0))
# Omegacor<- computeAssociations(m)
# supportLevel<- 0.95
# toPlot<- ((Omegacor[[1]]$support>supportLevel)+ (Omegacor[[1]]$support<(1-supportLevel))>0)*Omegacor[[1]]$mean
# corrplot(toPlot, method="color", col=colorRampPalette(c("blue", "white", "red"))(200))
par(mfrow=c(2,3),oma = c(1, 1, 1, 1))
corrplot(cor(hm_mod$Y), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Correlation cor(Y)")
corrplot(postRMean, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("R")
corrplot(Toplot_R, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Plot only non zero value")
corrplot(Tau_n, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Partial correlation matrix")
corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("True interactions")
Toplot_R
}
Rho_sign<-hm_inter(hm_mod)
# h_metric_e5<-metrics_hmsc(Rho_sign)
# cat("Success rate ")
# h_metric_e5$success_env
cat(" ")
cat(sprintf("Success rate: %s\n", metrics_hmsc(Rho_sign)$success_env))
## Success rate: 1
me10 <- load_object("model-2019-04-10-08-26-20.rda")
#load("model-2019-04-09-19-02-16.rda")
summary(me10)
## Summary for model '/tmp/RtmpKix4lZ/file40e957831178'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 683.734 minutes at time 2019-04-09 21:02:35.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
jsdm_conv(me10)
## Maximum Rhat value for Beta: 1.5602
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.7953
me10$mcmc.info[1:7]
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
j_metric_e10<-metrics_jsdm(me10)
cat("Success rate ")
## Success rate
j_metric_e10$success_env
## [1] 0.9333333
data<-sim_data$EnvEvenSp10
data <- list(
Y = subset(data, select = -env),
X = cbind(1, scale(poly(data$env, 2))),
covx = cov(cbind(1, scale(poly(data$env, 2)))),
K = 3,
J = ncol(data) - 1,
n = nrow(data),
I = diag(ncol(data) - 1),
df = ncol(data)
)
##########################################################################################
#Tau<-solve
#Tau_n<-to_prec(me10$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))
plot_cor_jsdm<-function(mod,y,interact=diag(ncol(y))){
par(mfrow=c(2,4),oma = c(3, 1, 2, 1))
corrplot(cor(y), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Correlation cor(Y)")
corrplot(mod$mean$EnvRho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("EnvRho")
corrplot(mod$mean$EnvRho*(!mod$overlap0$EnvRho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("EnvRho signif")
corrplot(mod$mean$Rho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Rho")
corrplot(mod$mean$Rho*(!mod$overlap0$Rho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Rho signif")
corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("True interactions")
}
plot_cor_jsdm(me10,data$Y)
data<-sim_data$EnvEvenSp10
#Rho_sign<-fit_gjam(data,5000,500,"./gjam_models/gjam10env.rda")
Rho_sign<-load_gjam(data,name="./gjam_models/gjam10env.rda")
#gje10<-load_object("./gjam_models/gjam10env.rda")
g_metric_e10<-metrics_gjam(Rho_sign)
cat("Success rate ")
## Success rate
g_metric_e10$success_env
## [1] 0.2666667
#to check posterior density of s in Sigma
#gje10<-load_object("./gjam_models/gjam10env.rda")
#plot(density(gje10$chains$sgibbs[,4]))
data<-sim_data$EnvEvenSp10
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm10env.rda" )
hm_conv(hm_mod)
Rho_sign<-hm_inter(hm_mod)
h_metric_e10<-metrics_hmsc(Rho_sign)
# cat("Success rate ")
# h_metric_e10$success_env
# h_metric_e10$success_fac
cat(" ")
cat(sprintf("Success rate: %s\n", h_metric_e10$success_env))
## Success rate: 1
cat(sprintf("Success rate: %s\n", h_metric_e10$success_fac))
me20 <- load_object("model-2019-04-11-19-06-02.rda")
#load("model-2019-04-09-19-02-16.rda")
summary(me20)
## Summary for model '/tmp/RtmpKix4lZ/file40e966e51ba7'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 2079.661 minutes at time 2019-04-10 08:26:21.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
me20$mcmc.info[1:7]
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
jsdm_conv(me20)
## Maximum Rhat value for Beta: 1.2174
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.3236
j_metric_e20<-metrics_jsdm(me20)
cat("Success rate ")
## Success rate
j_metric_e20$success_env
## [1] 0.9894737
data<-sim_data$EnvEvenSp20
data <- list(
Y = subset(data, select = -env),
X = cbind(1, scale(poly(data$env, 2))),
covx = cov(cbind(1, scale(poly(data$env, 2)))),
K = 3,
J = ncol(data) - 1,
n = nrow(data),
I = diag(ncol(data) - 1),
df = ncol(data)
)
##########################################################################################
plot_cor_jsdm(me20,data$Y)
data<-sim_data$EnvEvenSp20
#Rho_sign<-fit_gjam(data,5000,500,"./gjam_models/gjam20env.rda")
Rho_sign<-load_gjam(data,name="./gjam_models/gjam20env.rda")
#gje20<-load_object("./gjam_models/gjam20env.rda")
# g_metric_e20<-metrics_gjam(Rho_sign)
# cat("Success rate ")
# g_metric_e20$success_env
cat(" ")
cat(sprintf("Success rate env: %s\n", round(metrics_gjam(Rho_sign)$success_env,3)))
## Success rate env: 0.384
#to check posterior density of s in Sigma
#gje20<-load_object("./gjam_models/gjam20env.rda")
#plot(density(gje20$chains$sgibbs[,4]))
data<-sim_data$EnvEvenSp20
#fit_gjam(data,5000,500,"./gjam_models/gjam20env.rda")
#load_gjam(data,name="./gjam_models/gjam20env.rda")
#gje20<-load_object("./gjam_models/gjam20env.rda")
#to check posterior density of s in Sigma
#gje20<-load_object("./gjam_models/gjam20env.rda")
#plot(density(gje20$chains$sgibbs[,4]))
data <- list(
Y = subset(data, select = -env),
X = cbind(1, scale(poly(data$env, 2))),
covx = cov(cbind(1, scale(poly(data$env, 2)))),
K = 3,
J = ncol(data) - 1,
n = nrow(data),
I = diag(ncol(data) - 1),
df = ncol(data)
)
xdata<-as.data.frame(data$X[,-1])
colnames(xdata)<- c("env","env2")
ydata<-as.data.frame(data$Y)
###formula
rl <- list(r = 8, N = 20)
formula<-as.formula( ~env+ env2)
ml <- list(ng = 2500, burnin = 500, typeNames = 'PA', reductList = rl)
####fit
mod_gjam1 <- gjam(formula, xdata = xdata, ydata = ydata, modelList = ml)
##
## Observations and responses:
## [1] 500 20
## Warning in .setupReduct(modelList, S, Q, n): dimension reduction requires
## reductList$N < no. responses
##
## Dimension reduced from 20 X 20 -> 20 X 8 responses
## ===========================================================================
## expanding covariance chains
## ===========================================================================
save(mod_gjam1, file = "./gjam_models/gjam20env_dr.rda")
summary(mod_gjam1)
##
## Sensitivity by predictor variables f:
## Estimate SE CI_025 CI_975
## env 0.958 0.0836 0.823 1.16
## env2 1.310 0.0685 1.180 1.45
##
## Coefficient matrix B:
## sp01 sp02 sp03 sp04 sp05 sp06 sp07 sp08 sp09
## intercept -0.341 -0.101 0.0155 0.114 0.292 0.337 0.470 0.482 0.484
## env -0.473 -0.480 -0.4800 -0.474 -0.470 -0.460 -0.437 -0.383 -0.286
## env2 0.285 0.106 0.0122 -0.104 -0.255 -0.325 -0.438 -0.455 -0.478
## RMSPE 0.406 0.405 0.4090 0.411 0.404 0.403 0.409 0.423 0.434
## sp10 sp11 sp12 sp13 sp14 sp15 sp16 sp17
## intercept 0.4880 0.488 0.486 0.478 0.476 0.395 0.272 0.175
## env -0.0399 0.207 0.457 0.437 0.453 0.466 0.477 0.478
## env2 -0.4830 -0.474 -0.481 -0.478 -0.457 -0.354 -0.262 -0.175
## RMSPE 0.4430 0.446 0.428 0.425 0.405 0.407 0.407 0.409
## sp18 sp19 sp20
## intercept -0.00278 -0.092 -0.276
## env 0.48000 0.479 0.470
## env2 0.00499 0.092 0.237
## RMSPE 0.40400 0.407 0.404
##
## Coefficient matrix B:
## Estimate SE CI_025 CI_975 sig95
## sp01_intercept -0.34100 0.0696 -0.47300 -0.2060 *
## sp01_env -0.47300 0.0268 -0.49900 -0.3980 *
## sp01_env2 0.28500 0.0683 0.14900 0.4210 *
## sp02_intercept -0.10100 0.0722 -0.23800 0.0340
## sp02_env -0.48000 0.0235 -0.50000 -0.4170 *
## sp02_env2 0.10600 0.0674 -0.02310 0.2410
## sp03_intercept 0.01550 0.0734 -0.12500 0.1600
## sp03_env -0.48000 0.0210 -0.49900 -0.4240 *
## sp03_env2 0.01220 0.0677 -0.11900 0.1410
## sp04_intercept 0.11400 0.0751 -0.02790 0.2590
## sp04_env -0.47400 0.0306 -0.49900 -0.3990 *
## sp04_env2 -0.10400 0.0677 -0.23200 0.0286
## sp05_intercept 0.29200 0.0767 0.14100 0.4440 *
## sp05_env -0.47000 0.0288 -0.49900 -0.3930 *
## sp05_env2 -0.25500 0.0710 -0.39300 -0.1130 *
## sp06_intercept 0.33700 0.0692 0.19200 0.4660 *
## sp06_env -0.46000 0.0387 -0.49900 -0.3620 *
## sp06_env2 -0.32500 0.0692 -0.46300 -0.1900 *
## sp07_intercept 0.47000 0.0270 0.40000 0.4990 *
## sp07_env -0.43700 0.0667 -0.49900 -0.2660 *
## sp07_env2 -0.43800 0.0483 -0.49700 -0.3220 *
## sp08_intercept 0.48200 0.0176 0.43300 0.5000 *
## sp08_env -0.38300 0.0949 -0.49500 -0.1350 *
## sp08_env2 -0.45500 0.0369 -0.49900 -0.3630 *
## sp09_intercept 0.48400 0.0155 0.44200 0.5000 *
## sp09_env -0.28600 0.1140 -0.48700 -0.0525 *
## sp09_env2 -0.47800 0.0208 -0.49900 -0.4220 *
## sp10_intercept 0.48800 0.0116 0.45800 0.5000 *
## sp10_env -0.03990 0.1210 -0.25200 0.2070
## sp10_env2 -0.48300 0.0165 -0.50000 -0.4370 *
## sp11_intercept 0.48800 0.0116 0.45600 0.5000 *
## sp11_env 0.20700 0.1130 0.00925 0.4570 *
## sp11_env2 -0.47400 0.0241 -0.49900 -0.4070 *
## sp12_intercept 0.48600 0.0139 0.44800 0.5000 *
## sp12_env 0.45700 0.0400 0.34600 0.4990 *
## sp12_env2 -0.48100 0.0178 -0.50000 -0.4330 *
## sp13_intercept 0.47800 0.0210 0.42100 0.4990 *
## sp13_env 0.43700 0.0596 0.28400 0.4980 *
## sp13_env2 -0.47800 0.0214 -0.49900 -0.4210 *
## sp14_intercept 0.47600 0.0238 0.41300 0.4990 *
## sp14_env 0.45300 0.0529 0.31300 0.4990 *
## sp14_env2 -0.45700 0.0356 -0.49800 -0.3710 *
## sp15_intercept 0.39500 0.0608 0.26300 0.4920 *
## sp15_env 0.46600 0.0318 0.38100 0.4990 *
## sp15_env2 -0.35400 0.0682 -0.47700 -0.2120 *
## sp16_intercept 0.27200 0.0747 0.12700 0.4160 *
## sp16_env 0.47700 0.0259 0.40900 0.4990 *
## sp16_env2 -0.26200 0.0734 -0.40900 -0.1180 *
## sp17_intercept 0.17500 0.0756 0.02290 0.3210 *
## sp17_env 0.47800 0.0228 0.41300 0.4990 *
## sp17_env2 -0.17500 0.0720 -0.31000 -0.0364 *
## sp18_intercept -0.00278 0.0748 -0.14300 0.1420
## sp18_env 0.48000 0.0226 0.42700 0.4990 *
## sp18_env2 0.00499 0.0703 -0.13200 0.1430
## sp19_intercept -0.09200 0.0729 -0.23500 0.0550
## sp19_env 0.47900 0.0231 0.42000 0.4990 *
## sp19_env2 0.09200 0.0693 -0.04010 0.2280
## sp20_intercept -0.27600 0.0756 -0.42900 -0.1320 *
## sp20_env 0.47000 0.0320 0.39100 0.4990 *
## sp20_env2 0.23700 0.0698 0.10500 0.3780 *
##
## Last column indicates if 95% posterior distribution contains zero.
##
## Coefficient matrix B, standardized for X:
## NULL
##
## Last column indicates if 95% posterior distribution contains zero.
##
## Coefficient matrix B, standardized for X and W:
## NULL
##
## Last column indicates if 95% posterior distribution contains zero.
##
## Design Table
## env env2
## VIF 1 1
## env2 0 NA
##
## Sample contains n = 500 observations on S = 20 response variables, and 2 predictors. Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.415, and the DIC is 77360. Computation involved 2500 Gibbs steps, with a burnin of 500. Dimension reduction was implemented with N = 11 and r = 8.
Tau <- solve(mod_gjam1$parameters$sigMu)
Tau_n = to_prec(Tau)
#postH<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.95)
#postL<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.05)
#pH<-convert_to_m(postH)
#pL<-convert_to_m(postL)
#R_sign<-cov2cor(mod_gjam1$parameters$sigMu)*(!(pH>0 & pL<0))
g_metric_e20_dim_red<-metrics_gjam(Rho_sign)
cat("Success rate ")
## Success rate
g_metric_e20_dim_red$success_env
## [1] 0.3842105
par(mfrow=c(2,3),oma = c(1, 1, 1, 1))
corrplot(cor(ydata), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Correlation cor(Y)")
corrplot(mod_gjam1$parameters$corMu, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("R")
corrplot(mod_gjam1$parameters$ematrix, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("E matrix")
# corrplot(Tau_n, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
# title("Tau")
# corrplot(R_sign, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
# title("R signif")
corrplot(diag(20), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("True interactions")
data<-sim_data$EnvEvenSp20
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm20env.rda" )
hm_conv(hm_mod)
Rho_sign<-hm_inter(hm_mod)
# h_metric_e20<-metrics_hmsc(Rho_sign)
# cat("Success score ")
# h_metric_e20$success_env
cat(" ")
cat(sprintf("Success rate non-interacting: %s\n", round(metrics_hmsc(Rho_sign)$success_env,3)))
## Success rate non-interacting: 1
mfd5 <- load_object("model-2019-04-11-19-35-11.rda")
summary(mfd5)
## Summary for model '/tmp/RtmpKix4lZ/file40e94e482135'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 29.13 minutes at time 2019-04-11 19:06:03.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
#mfd5$Rhat
jsdm_conv(mfd5)
## Maximum Rhat value for Beta: 1.3081
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.1421
mfd5$mcmc.info[1:7]
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
j_metric_facDense5<-metrics_jsdm(mfd5,fac = fac_inter[[4]],only_env = F)
cat("Success rate ")
## Success rate
j_metric_facDense5$success_env
## [1] 1
j_metric_facDense5$success_fac
## [1] 0
data<-sim_data$FacDenseSp5
data <- list(
Y = subset(data, select = -env),
X = cbind(1, scale(poly(data$env, 2))),
covx = cov(cbind(1, scale(poly(data$env, 2)))),
K = 3,
J = ncol(data) - 1,
n = nrow(data),
I = diag(ncol(data) - 1),
df = ncol(data)
)
##########################################################################################
#Tau<-solve
#Tau_n<-to_prec(me10$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))
plot_cor_jsdm(mfd5,data$Y,fac_inter[[4]])
data<-sim_data$FacDenseSp5
#Rho_sign<-fit_gjam(data,10000,2000,"./gjam_models/gjam5f.rda",interact=fac_inter[[4]])
Rho_sign<-load_gjam(data,name="./gjam_models/gjam5f.rda", interact=fac_inter[[4]])
#gjfd5<-load_object("./gjam_models/gjam5f.rda")
g_metric_facDense5<-metrics_gjam(Rho_sign,fac=fac_inter[[4]], only_env = F)
#cat("Success rate ")
#g_metric_facDense5$success_env
#g_metric_facDense5$success_fac
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam5f.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facDense5$success_env,3)))
## Success rate for non-interacting: 0.375
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facDense5$success_fac,3)))
## Success rate for facilitation: 0.5
data<-sim_data$FacDenseSp5
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm5fd.rda" )
hm_conv(hm_mod)
hm_inter(hm_mod, interact = fac_inter[[4]])
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 0 0 0
## [2,] 0 1 0 0 0
## [3,] 0 0 1 0 0
## [4,] 0 0 0 1 0
## [5,] 0 0 0 0 1
Rho_sign<-hm_inter(hm_mod)
h_metric_facDense5<-metrics_hmsc(Rho_sign,fac=fac_inter[[4]],only_env = F)
#cat("Success rate ")
# h_metric_facDense5$success_env
# h_metric_facDense5$success_fac
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facDense5$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_facDense5$success_fac,3)))
## Success rate for facilitation: 0
## Summary for model '/tmp/RtmpKix4lZ/file40e96d09e31e'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 684.487 minutes at time 2019-04-11 19:35:12.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.2866
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4808
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate
## [1] 0.925
## [1] 0
data<-sim_data$FacDenseSp10
#Rho_sign<-fit_gjam(data,5000,500,"./gjam_models/gjam10fd.rda",interact=fac_inter[[5]])
load_gjam(data,name="./gjam_models/gjam10fd.rda", interact=fac_inter[[5]])
## sp01 sp02 sp03 sp04 sp05 sp06
## sp01 1.0000000 -0.3904084 0.0000000 0.2493194 -0.1879177 -0.6181350
## sp02 -0.3904084 1.0000000 0.0000000 -0.3751301 0.0000000 0.0000000
## sp03 0.0000000 0.0000000 1.0000000 0.2974028 0.0000000 0.0000000
## sp04 0.2493194 -0.3751301 0.2974028 1.0000000 0.0000000 0.3768348
## sp05 -0.1879177 0.0000000 0.0000000 0.0000000 1.0000000 0.0000000
## sp06 -0.6181350 0.0000000 0.0000000 0.3768348 0.0000000 1.0000000
## sp07 0.0000000 0.0000000 0.0000000 -0.7225666 -0.1509510 -0.3458261
## sp08 0.0000000 0.4069187 0.0000000 0.0000000 0.0000000 0.0000000
## sp09 -0.1450869 0.0000000 -0.1989728 0.0000000 0.0000000 0.0000000
## sp10 -0.5235049 -0.2648058 0.0000000 0.0000000 0.0000000 0.6694006
## sp07 sp08 sp09 sp10
## sp01 0.0000000 0.0000000 -0.1450869 -0.5235049
## sp02 0.0000000 0.4069187 0.0000000 -0.2648058
## sp03 0.0000000 0.0000000 -0.1989728 0.0000000
## sp04 -0.7225666 0.0000000 0.0000000 0.0000000
## sp05 -0.1509510 0.0000000 0.0000000 0.0000000
## sp06 -0.3458261 0.0000000 0.0000000 0.6694006
## sp07 1.0000000 0.0000000 0.3544135 0.3680399
## sp08 0.0000000 1.0000000 0.0000000 0.0000000
## sp09 0.3544135 0.0000000 1.0000000 0.3121062
## sp10 0.3680399 0.0000000 0.3121062 1.0000000
#gjfd5<-load_object("./gjam_models/gjam10fd.rda")
g_metric_facDense10<-metrics_gjam(Rho_sign,fac=fac_inter[[5]], only_env = F)
# cat("Success rate ")
# g_metric_facDense10$success_env
# g_metric_facDense10$success_fac
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam10fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facDense10$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facDense10$success_fac,3)))
## Success rate for facilitation: 0
## Success rate for non-interacting: 1
## Success rate for facilitation: 0
## Summary for model '/tmp/RtmpKix4lZ/file40e91ec9ff9'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 2062.467 minutes at time 2019-04-12 06:59:42.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.6132
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.508
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate
## [1] 0.9833333
## [1] 0
data<-sim_data$FacDenseSp20
#Rho_sign<-fit_gjam(data,10000,1500,"./gjam_models/gjam20fd.rda",interact=fac_inter[[6]])
Rho_sign<-load_gjam(data,name="./gjam_models/gjam20fd.rda", interact=fac_inter[[6]])
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
g_metric_facDense20<-metrics_gjam(Rho_sign,fac=fac_inter[[6]], only_env = F)
# cat("Success rate ")
# g_metric_facDense20$success_env
# g_metric_facDense20$success_fac
#
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facDense20$success_env,3)))
## Success rate for non-interacting: 0.411
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facDense20$success_fac,3)))
## Success rate for facilitation: 0.5
data<-sim_data$FacDenseSp20
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm20fd.rda" )
hm_conv(hm_mod)
Rho_sign<-hm_inter(hm_mod, interact = fac_inter[[6]])
h_metric_facDense20<-metrics_hmsc(Rho_sign,fac=fac_inter[[6]],only_env = F)
# cat("Success rate ")
# h_metric_facDense20$success_env
# h_metric_facDense20$success_fac
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facDense20$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_facDense20$success_fac,3)))
## Success rate for facilitation: 0
## Summary for model '/tmp/RtmpKix4lZ/file40e917d13dae'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 29.05 minutes at time 2019-04-13 17:22:13.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.373
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4048
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.889
## Success rate for facilitation: 1
data<-sim_data$FacSparseSp5
#Rho_sign<-fit_gjam(data,2500,500,"./gjam_models/gjam5fs.rda",interact=fac_inter[[7]])
Rho_sign<-load_gjam(data,name="./gjam_models/gjam5fs.rda", interact=fac_inter[[7]])
#gjfd5<-load_object("./gjam_models/gjam5fs.rda")
g_metric_facSparse5<-metrics_gjam(Rho_sign,fac=fac_inter[[7]], only_env = F)
# cat("Success rate ")
# g_metric_facSparse5$success_env
# g_metric_facSparse5$success_fac
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facSparse5$success_env,3)))
## Success rate for non-interacting: 0.556
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facSparse5$success_fac,3)))
## Success rate for facilitation: 0
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
data<-sim_data$FacSparseSp5
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm5fs.rda" )
hm_conv(hm_mod)
Rho_sign<-hm_inter(hm_mod, interact = fac_inter[[7]])
h_metric_facSparse5<-metrics_hmsc(Rho_sign,fac=fac_inter[[8]],only_env = F)
# cat("Success rate ")
# h_metric_facSparse5$success_env
# h_metric_facSparse5$success_fac
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facSparse5$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_facSparse5$success_fac,3)))
## Success rate for facilitation: 0
## Summary for model '/tmp/RtmpKix4lZ/file40e9113ddc9'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 686.687 minutes at time 2019-04-13 17:51:16.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.5479
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4746
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.93
## Success rate for facilitation: 0
data<-sim_data$FacSparseSp10
#Rho_sign<-fit_gjam(data,2500,500,"./gjam_models/gjam10fs.rda",interact=fac_inter[[8]])
Rho_sign<-load_gjam(data,name="./gjam_models/gjam10fs.rda", interact=fac_inter[[8]])
#gjfd5<-load_object("./gjam_models/gjam10fs.rda")
g_metric_facSparse10<-metrics_gjam(Rho_sign,fac=fac_inter[[8]], only_env = F)
# cat("Success rate ")
# g_metric_facSparse10$success_env
# g_metric_facSparse10$success_fac
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam10fs.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facSparse10$success_env,3)))
## Success rate for non-interacting: 0.372
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facSparse10$success_fac,3)))
## Success rate for facilitation: 0
data<-sim_data$FacSparseSp10
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm10fs.rda" )
hm_conv(hm_mod)
Rho_sign<-hm_inter(hm_mod, interact = fac_inter[[8]])
h_metric_facSparse10<-metrics_hmsc(Rho_sign,fac=fac_inter[[8]],only_env = F)
# cat("Success rate ")
# h_metric_facSparse10$success_env
# h_metric_facSparse10$success_fac
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facSparse10$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_facSparse10$success_fac,3)))
## Success rate for facilitation: 0
## Summary for model '/tmp/RtmpKix4lZ/file40e9629c97ad'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 2056.316 minutes at time 2019-04-14 05:17:59.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.5653
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4826
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.978
## Success rate for facilitation: 0
data<-sim_data$FacSparseSp20
#fit_gjam(data,10000,1500,"./gjam_models/gjam20fs.rda",interact=fac_inter[[9]])
load_gjam(data,name="./gjam_models/gjam20fs.rda", interact=fac_inter[[9]])
## sp01 sp02 sp03 sp04 sp05 sp06
## sp01 1.0000000 0.0000000 0.0000000 0.2481647 -0.1405985 -0.1778950
## sp02 0.0000000 1.0000000 -0.1750012 0.0000000 0.0000000 -0.6945088
## sp03 0.0000000 -0.1750012 1.0000000 0.0000000 0.1633683 0.3398000
## sp04 0.2481647 0.0000000 0.0000000 1.0000000 0.1812365 0.0000000
## sp05 -0.1405985 0.0000000 0.1633683 0.1812365 1.0000000 0.0000000
## sp06 -0.1778950 -0.6945088 0.3398000 0.0000000 0.0000000 1.0000000
## sp07 -0.2262923 0.0000000 0.2351144 -0.2474855 -0.3272462 0.0000000
## sp08 0.0000000 0.4861375 0.0000000 -0.2295886 0.0000000 -0.3307800
## sp09 0.0000000 0.0000000 0.5327952 0.0000000 0.0000000 0.3040397
## sp10 0.2989764 0.0000000 0.4502637 0.0000000 0.0000000 -0.1168824
## sp11 0.3747704 -0.2002289 0.3865259 0.2438452 -0.2970431 0.0000000
## sp12 -0.2809890 0.0000000 0.1516279 -0.4507748 0.0000000 -0.1320808
## sp13 0.0000000 -0.2293087 0.3528826 0.0000000 0.2622547 0.2939977
## sp14 0.0000000 0.1499564 0.0000000 0.0000000 0.0000000 0.0000000
## sp15 0.3873898 0.0000000 -0.2520398 -0.2776349 0.1619795 0.0000000
## sp16 -0.2582433 0.0000000 -0.1597638 -0.2483336 -0.2604524 0.2649414
## sp17 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.2525415
## sp18 0.0000000 -0.1554376 0.0000000 0.0000000 0.0000000 0.0000000
## sp19 -0.1567031 0.2745205 0.1115394 0.0000000 0.0000000 0.0000000
## sp20 -0.3296659 0.0000000 0.0000000 0.0000000 0.1530916 -0.2474593
## sp07 sp08 sp09 sp10 sp11 sp12
## sp01 -0.2262923 0.0000000 0.0000000 0.2989764 0.3747704 -0.2809890
## sp02 0.0000000 0.4861375 0.0000000 0.0000000 -0.2002289 0.0000000
## sp03 0.2351144 0.0000000 0.5327952 0.4502637 0.3865259 0.1516279
## sp04 -0.2474855 -0.2295886 0.0000000 0.0000000 0.2438452 -0.4507748
## sp05 -0.3272462 0.0000000 0.0000000 0.0000000 -0.2970431 0.0000000
## sp06 0.0000000 -0.3307800 0.3040397 -0.1168824 0.0000000 -0.1320808
## sp07 1.0000000 0.1866792 0.1815724 -0.1164417 0.0000000 -0.1966267
## sp08 0.1866792 1.0000000 0.2297215 0.0000000 0.0000000 0.0000000
## sp09 0.1815724 0.2297215 1.0000000 0.1862230 0.0000000 0.0000000
## sp10 -0.1164417 0.0000000 0.1862230 1.0000000 0.1139547 0.0000000
## sp11 0.0000000 0.0000000 0.0000000 0.1139547 1.0000000 0.0000000
## sp12 -0.1966267 0.0000000 0.0000000 0.0000000 0.0000000 1.0000000
## sp13 0.0000000 -0.4251175 -0.1844692 -0.1109607 0.0000000 0.0000000
## sp14 0.0000000 -0.3169754 0.0000000 -0.1806383 0.1229519 0.0000000
## sp15 -0.2393655 0.3940648 0.0000000 0.0000000 -0.1572325 0.0000000
## sp16 0.0000000 0.0000000 -0.2478583 -0.1822943 0.0000000 0.0000000
## sp17 -0.4179798 -0.2271328 0.2818370 -0.3565917 -0.2899321 0.0000000
## sp18 -0.3452972 0.0000000 0.3590445 -0.1968948 0.0000000 0.0000000
## sp19 0.0000000 -0.3549123 -0.3927233 0.0000000 0.0000000 0.2486426
## sp20 -0.2687936 0.0000000 0.0000000 -0.1423872 0.3274940 0.4736903
## sp13 sp14 sp15 sp16 sp17 sp18
## sp01 0.0000000 0.0000000 0.3873898 -0.2582433 0.0000000 0.0000000
## sp02 -0.2293087 0.1499564 0.0000000 0.0000000 0.0000000 -0.1554376
## sp03 0.3528826 0.0000000 -0.2520398 -0.1597638 0.0000000 0.0000000
## sp04 0.0000000 0.0000000 -0.2776349 -0.2483336 0.0000000 0.0000000
## sp05 0.2622547 0.0000000 0.1619795 -0.2604524 0.0000000 0.0000000
## sp06 0.2939977 0.0000000 0.0000000 0.2649414 0.2525415 0.0000000
## sp07 0.0000000 0.0000000 -0.2393655 0.0000000 -0.4179798 -0.3452972
## sp08 -0.4251175 -0.3169754 0.3940648 0.0000000 -0.2271328 0.0000000
## sp09 -0.1844692 0.0000000 0.0000000 -0.2478583 0.2818370 0.3590445
## sp10 -0.1109607 -0.1806383 0.0000000 -0.1822943 -0.3565917 -0.1968948
## sp11 0.0000000 0.1229519 -0.1572325 0.0000000 -0.2899321 0.0000000
## sp12 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## sp13 1.0000000 0.1641325 0.0000000 0.0000000 0.0000000 -0.1065072
## sp14 0.1641325 1.0000000 0.0000000 0.2800982 0.2165972 0.0000000
## sp15 0.0000000 0.0000000 1.0000000 0.0000000 0.2232036 0.0000000
## sp16 0.0000000 0.2800982 0.0000000 1.0000000 0.0000000 0.2704149
## sp17 0.0000000 0.2165972 0.2232036 0.0000000 1.0000000 0.2673245
## sp18 -0.1065072 0.0000000 0.0000000 0.2704149 0.2673245 1.0000000
## sp19 0.6487585 0.0000000 -0.2207364 0.0000000 0.0000000 -0.3499635
## sp20 -0.2588774 -0.4217491 -0.2976427 0.0000000 0.0000000 0.2694292
## sp19 sp20
## sp01 -0.1567031 -0.3296659
## sp02 0.2745205 0.0000000
## sp03 0.1115394 0.0000000
## sp04 0.0000000 0.0000000
## sp05 0.0000000 0.1530916
## sp06 0.0000000 -0.2474593
## sp07 0.0000000 -0.2687936
## sp08 -0.3549123 0.0000000
## sp09 -0.3927233 0.0000000
## sp10 0.0000000 -0.1423872
## sp11 0.0000000 0.3274940
## sp12 0.2486426 0.4736903
## sp13 0.6487585 -0.2588774
## sp14 0.0000000 -0.4217491
## sp15 -0.2207364 -0.2976427
## sp16 0.0000000 0.0000000
## sp17 0.0000000 0.0000000
## sp18 -0.3499635 0.2694292
## sp19 1.0000000 0.0000000
## sp20 0.0000000 1.0000000
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
g_metric_facSparse20<-metrics_gjam(Rho_sign,fac=fac_inter[[9]], only_env = F)
# cat("Success rate ")
# g_metric_facSparse20$success_env
# g_metric_facSparse20$success_fac
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facSparse20$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facSparse20$success_fac,3)))
## Success rate for facilitation: 0
data<-sim_data$FacSparseSp20
hm_mod<-fit_hmsc(data,nsamples=3000, nchains=2,"Load",name="./HMmodels/hm20fs.rda" )
hm_conv(hm_mod)
hm_inter(hm_mod, nsamples=3000, nchains=2,interact = fac_inter[[9]])
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 1 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 1 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 1 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 1 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 1 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 1 0 0 0 0 0 0 0
## [7,] 0 0 0 0 0 0 1 0 0 0 0 0 0
## [8,] 0 0 0 0 0 0 0 1 0 0 0 0 0
## [9,] 0 0 0 0 0 0 0 0 1 0 0 0 0
## [10,] 0 0 0 0 0 0 0 0 0 1 0 0 0
## [11,] 0 0 0 0 0 0 0 0 0 0 1 0 0
## [12,] 0 0 0 0 0 0 0 0 0 0 0 1 0
## [13,] 0 0 0 0 0 0 0 0 0 0 0 0 1
## [14,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [15,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [16,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [17,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [18,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [19,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [20,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## [1,] 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0
## [7,] 0 0 0 0 0 0 0
## [8,] 0 0 0 0 0 0 0
## [9,] 0 0 0 0 0 0 0
## [10,] 0 0 0 0 0 0 0
## [11,] 0 0 0 0 0 0 0
## [12,] 0 0 0 0 0 0 0
## [13,] 0 0 0 0 0 0 0
## [14,] 1 0 0 0 0 0 0
## [15,] 0 1 0 0 0 0 0
## [16,] 0 0 1 0 0 0 0
## [17,] 0 0 0 1 0 0 0
## [18,] 0 0 0 0 1 0 0
## [19,] 0 0 0 0 0 1 0
## [20,] 0 0 0 0 0 0 1
h_metric_facSparse20<-metrics_hmsc(Rho_sign,fac=fac_inter[[9]],only_env = F)
# cat("Success rate ")
# h_metric_facSparse20$success_env
# h_metric_facSparse20$success_fac
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facSparse20$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_facSparse20$success_fac,3)))
## Success rate for facilitation: 0
## Summary for model '/tmp/RtmpKix4lZ/file40e9127bbb96'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 29.308 minutes at time 2019-04-15 15:34:20.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## Successful convergence based on Rhat values (all < 1.1).
## Maximum Rhat value for Beta: 1.0324
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.0313
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate
## [1] 1
## [1] 1
## Success rate for non-interacting: 1
## Success rate for facilitation: 0
data<-sim_data$CompDenseSp5
#Rho_sign<-fit_gjam(data,10000,1000,"./gjam_models/gjam5cmpd.rda",interact=comp_inter[[10]])
Rho_sign<-load_gjam(data,name="./gjam_models/gjam5cmpd.rda", interact=(-1)*comp_inter[[10]])
#gjfd5<-load_object("./gjam_models/gjam5cmpd.rda")
g_metric_compDense5<-metrics_gjam(Rho_sign,comp=comp_inter[[10]], only_env = F)
# cat("Success rate ")
# g_metric_compDense5$success_env
# g_metric_compDense5$success_comp
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compDense5$success_env,3)))
## Success rate for non-interacting: 0.625
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_compDense5$success_comp,3)))
## Success rate for facilitation: 1
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam5cmpd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
data<-sim_data$CompDenseSp5
hm_mod<-fit_hmsc(data,"Load",nsamples=3000, nchains=2,name="./HMmodels/hm5cmpd.rda" )
#hm_mod<-fit_hmsc(data,"Fit",nsamples=3000, nchains=2 )
hm_conv(hm_mod)
Rho_sign<-hm_inter(hm_mod, nsamples=3000, nchains=2,interact = (-1)*comp_inter[[10]])
h_metric_compDense5<-metrics_hmsc(Rho_sign,comp=comp_inter[[10]],only_env = F)
# cat("Success rate ")
# h_metric_compDense5$success_env
# h_metric_compDense5$success_comp
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compDense5$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_compDense5$success_comp,3)))
## Success rate for facilitation: 1
## Summary for model '/tmp/RtmpKix4lZ/file40e9559ce9d8'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 695.614 minutes at time 2019-04-15 16:03:39.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.1845
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.1771
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.875
## Success rate for facilitation: 1
data<-sim_data$CompDenseSp10
#Rho_sign<-fit_gjam(data,5000,500,"./gjam_models/gjam10cmpd.rda",interact=comp_inter[[11]])
Rho_sign<-load_gjam(data,name="./gjam_models/gjam10cmpd.rda", interact= (-1)*comp_inter[[11]])
#gjfd5<-load_object("./gjam_models/gjam10cmpd.rda")
g_metric_compDense10<-metrics_gjam(Rho_sign,comp=comp_inter[[11]], only_env = F)
# cat("Success rate ")
# g_metric_compDense10$success_env
# g_metric_compDense10$success_comp
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compDense10$success_env,3)))
## Success rate for non-interacting: 0.35
cat(sprintf("Success rate for competition: %s\n", round(g_metric_compDense10$success_comp,3)))
## Success rate for competition: 1
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
data<-sim_data$CompDenseSp10
hm_mod<-fit_hmsc(data,"Load",nsamples=5000, nchains=2,name="./HMmodels/hm10cmpd.rda" )
hm_conv(hm_mod)
Rho_sign<-hm_inter(hm_mod,nsamples=5000, nchains=2, interact = (-1)*comp_inter[[11]])
h_metric_compDense10<-metrics_hmsc(Rho_sign,comp=comp_inter[[11]],only_env = F)
# cat("Success rate ")
# h_metric_compDense10$success_env
# h_metric_compDense10$success_comp
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compDense10$success_env,3)))
## Success rate for non-interacting: 0.875
cat(sprintf("Success rate for competition: %s\n", round(h_metric_compDense10$success_comp,3)))
## Success rate for competition: 1
## Summary for model '/tmp/RtmpKix4lZ/file40e96843124a'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 2061.424 minutes at time 2019-04-16 03:39:17.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.9065
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.6944
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.956
## Success rate for competition: 0.7
data<-sim_data$CompDenseSp20
#Rho_sign<-fit_gjam(data,10000,1000,"./gjam_models/gjam20cmpd.rda",interact= (-1)*comp_inter[[12]])
Rho_sign<-load_gjam(data,name="./gjam_models/gjam20cmpd.rda", interact= (-1)*comp_inter[[12]])
#gjfd5<-load_object("./gjam_models/gjam20cmpd.rda")
g_metric_compDense20<-metrics_gjam(Rho_sign,comp=comp_inter[[12]], only_env = F)
# cat("Success rate ")
# g_metric_compDense20$success_env
# g_metric_compDense20$success_comp
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compDense20$success_env,3)))
## Success rate for non-interacting: 0.406
cat(sprintf("Success rate for competition: %s\n", round(g_metric_compDense20$success_comp,3)))
## Success rate for competition: 0.9
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
data<-sim_data$CompDenseSp20
hm_mod<-fit_hmsc(data,"Load",nsamples=5000, nchains=2,name="./HMmodels/hm20cmpd.rda" )
hm_conv(hm_mod)
Rho_sign<-hm_inter(hm_mod, nsamples=5000, nchains=2,interact = (-1)*comp_inter[[12]])
h_metric_compDense20<-metrics_hmsc(Rho_sign,comp=comp_inter[[12]],only_env = F)
# cat("Success rate ")
# h_metric_compDense20$success_env
# h_metric_compDense20$success_comp
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compDense20$success_env,3)))
## Success rate for non-interacting: 0.967
cat(sprintf("Success rate for competition: %s\n", round(h_metric_compDense20$success_comp,3)))
## Success rate for competition: 0.4
## Summary for model '/tmp/RtmpKix4lZ/file40e97e225117'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 29.251 minutes at time 2019-04-17 14:00:46.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.1876
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.2038
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate for non-interacting: 1
## Success rate for competition: 1
data<-sim_data$CompSparseSp5
#Rho_sign<-fit_gjam(data,5000,2000,"./gjam_models/gjam5cmps.rda",interact= (-1)*comp_inter[[13]])
Rho_sign<-load_gjam(data,name="./gjam_models/gjam5cmps.rda", interact= (-1)*comp_inter[[13]])
#gjfd5<-load_object("./gjam_models/gjam5cmps.rda")
g_metric_compSparse5<-metrics_gjam(Rho_sign,comp=comp_inter[[13]], only_env = F)
cat("Success rate ")
## Success rate
g_metric_compSparse5$success_env
## [1] 0.7777778
g_metric_compSparse5$success_comp
## [1] 1
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compSparse5$success_env,3)))
## Success rate for non-interacting: 0.778
cat(sprintf("Success rate for competition: %s\n", round(g_metric_compSparse5$success_comp,3)))
## Success rate for competition: 1
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam5cmps.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
data<-sim_data$CompSparseSp5
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hm5cmps.rda" )
hm_conv(hm_mod)
Rho_sign<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact = (-1)*comp_inter[[13]])
h_metric_facSparse5<-metrics_hmsc(Rho_sign,comp=comp_inter[[13]],only_env = F)
# cat("Success rate ")
# h_metric_facSparse5$success_env
# h_metric_facSparse5$success_comp
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facSparse5$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(h_metric_facSparse5$success_comp,3)))
## Success rate for competition: 1
## Summary for model '/tmp/RtmpKix4lZ/file40e978d641ed'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 685.719 minutes at time 2019-04-17 14:30:01.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.116
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.1753
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate
## [1] 0.9534884
## [1] 1
## Success rate for non-interacting: 0.953
## Success rate for competition: 1
data<-sim_data$CompSparseSp10
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjam10cmps.rda",interact= (-1)*comp_inter[[14]])
Rho_sign<-load_gjam(data,name="./gjam_models/gjam10cmps.rda", interact= (-1)*comp_inter[[14]])
g_metric_compSparse10<-metrics_gjam(Rho_sign,comp=comp_inter[[14]], only_env = F)
# cat("Success rate ")
# g_metric_compSparse10$success_env
# g_metric_compSparse10$success_comp
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compSparse10$success_env,3)))
## Success rate for non-interacting: 0.349
cat(sprintf("Success rate for competition: %s\n", round(g_metric_compSparse10$success_comp,3)))
## Success rate for competition: 1
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam10cmps.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
data<-sim_data$CompSparseSp10
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hm10cmps.rda" )
hm_conv(hm_mod)
Rho_sign<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact = (-1)*comp_inter[[14]])
h_metric_facSparse10<-metrics_hmsc(Rho_sign,comp=comp_inter[[14]],only_env = F)
cat("Success rate ")
## Success rate
h_metric_facSparse10$success_env
## [1] 0.9534884
h_metric_facSparse10$success_comp
## [1] 0.5
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facSparse10$success_env,3)))
## Success rate for non-interacting: 0.953
cat(sprintf("Success rate for competition: %s\n", round(h_metric_facSparse10$success_comp,3)))
## Success rate for competition: 0.5
## Summary for model '/tmp/RtmpKix4lZ/file40e9379ee963'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 2084.707 minutes at time 2019-04-18 01:55:46.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.2652
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.3883
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.968
## Success rate for competition: 0
data<-sim_data$CompSparseSp20
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjam20cmps.rda",interact= (-1)*comp_inter[[15]])
Rho_sign<-load_gjam(data,name="./gjam_models/gjam20cmps.rda", interact= (-1)*comp_inter[[15]])
g_metric_compSparse20<-metrics_gjam(Rho_sign,comp=comp_inter[[15]], only_env = F)
# cat("Success rate ")
# g_metric_compSparse20$success_env
# g_metric_compSparse20$success_comp
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compSparse20$success_env,3)))
## Success rate for non-interacting: 0.296
cat(sprintf("Success rate for competition: %s\n", round(g_metric_compSparse20$success_comp,3)))
## Success rate for competition: 0.75
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam20cmps.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
data<-sim_data$CompSparseSp20
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hm20cmps.rda" )
hm_conv(hm_mod)
Rho_sign<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact = (-1)*comp_inter[[15]])
h_metric_compSparse20<-metrics_hmsc(Rho_sign,comp=comp_inter[[15]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compSparse20$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(h_metric_compSparse20$success_comp,3)))
## Success rate for competition: 0
## Summary for model '/tmp/RtmpKix4lZ/file40e93b49ad29'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 28.267 minutes at time 2019-04-19 12:40:31.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.1667
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.2086
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.5
## Success rate for facilitation: 0.5
data<-sim_data$FacCompDenseSp5
#Rho_sign<-fit_gjam(data,30000,15000,"./gjam_models/gjam5cmp_facd.rda",interact= (-1)*comp_inter[[16]]+fac_inter[[16]])
Rho_sign2<-load_gjam(data,"./gjam_models/gjam5cmp_facd2.rda",interact= (-1)*comp_inter[[16]]+fac_inter[[16]])
#Rho_sign3<-fit_gjam(data,5000,1000,"./gjam_models/gjam5cmp_facd2.rda",interact= (-1)*comp_inter[[16]]+fac_inter[[16]])
#Rho_sign4<-fit_gjam(data,5000,2500,"./gjam_models/gjam5cmp_facd.rda",interact= (-1)*comp_inter[[16]]+fac_inter[[16]])
#Rho_sign5<-fit_gjam(data,20000,5000,"./gjam_models/gjam5cmp_facd2.rda",interact= (-1)*comp_inter[[16]]+fac_inter[[16]])
#Rho_sign<-load_gjam(data,name="./gjam_models/gjam5cmp_facd.rda", (-1)*comp_inter[[16]]+fac_inter[[16]])
g_metric_cmp_facd5<-metrics_gjam(Rho_sign2,comp=comp_inter[[16]], fac=fac_inter[[16]], only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_cmp_facd5$success_env,3)))
## Success rate for non-interacting: 0.833
cat(sprintf("Success rate for competition: %s\n", round(g_metric_cmp_facd5$success_comp,3)))
## Success rate for competition: 1
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_cmp_facd5$success_fac,3)))
## Success rate for facilitation: 1
# par(mfrow=c(2,3))
#corrplot((Rho_sign+Rho_sign2)/2, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", col=cols(200),type = "lower")
# corrplot(Rho_sign2, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", col=cols(200),type = "lower")
# corrplot(Rho_sign3, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", col=cols(200),type = "lower")
# corrplot(Rho_sign4, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", col=cols(200),type = "lower")
# corrplot(Rho_sign5, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", col=cols(200),type = "lower")
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam20cmps.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
data<-sim_data$FacCompDenseSp5
hm_mod<-fit_hmsc(data,"Load",nsamples=3000, nchains=2,name="./HMmodels/hmcmp_facd5.rda" )
hm_conv(hm_mod)
Rho_sign<-hm_inter(hm_mod, nsamples=3000, nchains=2,interact = (-1)*comp_inter[[16]]+ fac_inter[[16]])
h_metric_FacCompDenseSp5<-metrics_hmsc(Rho_sign,comp=comp_inter[[16]],fac= fac_inter[[16]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompDenseSp5$success_env,3)))
## Success rate for non-interacting: 0.833
cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompDenseSp5$success_comp,3)))
## Success rate for competition: 0.5
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompDenseSp5$success_fac,3)))
## Success rate for facilitation: 0.5
data<-sim_data$FacCompDenseSp10
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjamcmp_facd10.rda",interact= (-1)*comp_inter[[17]]+fac_inter[[17]])
Rho_sign<-load_gjam(data,name="./gjam_models/gjamcmp_facd10.rda", interact= (-1)*comp_inter[[17]]+fac_inter[[17]])
g_metric_FacCompDenseSp10<-metrics_gjam(Rho_sign,comp=comp_inter[[17]], fac=fac_inter[[17]], only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_FacCompDenseSp10$success_env,3)))
## Success rate for non-interacting: 0.343
cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompDenseSp10$success_comp,3)))
## Success rate for competition: 0.6
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_FacCompDenseSp10$success_fac,3)))
## Success rate for facilitation: 0.667
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hmcmp_facd10.rda" )
hm_conv(hm_mod)
Rho_sign<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact = (-1)*comp_inter[[17]] + fac_inter[[17]])
h_metric_FacCompDenseSp10<-metrics_hmsc(Rho_sign,comp=comp_inter[[17]],fac=fac_inter[[17]], only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompDenseSp10$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompDenseSp10$success_comp,3)))
## Success rate for competition: 0
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompDenseSp10$success_fac,3)))
## Success rate for facilitation: 0
data<-sim_data$FacCompDenseSp20
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hmcmp_facd20.rda" )
hm_conv(hm_mod)
Rho_sign<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact = (-1)*comp_inter[[18]] + fac_inter[[18]])
h_metric_FacCompDenseSp20<-metrics_hmsc(Rho_sign,comp=comp_inter[[18]],fac=fac_inter[[18]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompDenseSp20$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompDenseSp20$success_comp,3)))
## Success rate for competition: 0
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompDenseSp20$success_fac,3)))
## Success rate for facilitation: 0
data<-sim_data$FacCompSparseSp5
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjamcmp_facs5.rda",interact= (-1)*comp_inter[[19]]+ fac_inter[[19]] )
Rho_sign<-load_gjam(data,name="./gjam_models/gjamcmp_facs5.rda", interact= (-1)*comp_inter[[19]]+ fac_inter[[19]])
g_metric_FacCompSparseSp5<-metrics_gjam(Rho_sign,comp=comp_inter[[15]], only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_FacCompSparseSp5$success_env,3)))
## Success rate for non-interacting: 0.778
cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompSparseSp5$success_comp,3)))
## Success rate for competition: 1
cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompSparseSp5$success_comp,3)))
## Success rate for competition: 1
data<-sim_data$FacCompSparseSp5
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hmcmp_facs5.rda" )
hm_conv(hm_mod)
Rho_sign<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact = (-1)*comp_inter[[19]]+ fac_inter[[19]])
FacCompSparseSp5<-metrics_hmsc(Rho_sign,comp=comp_inter[[19]],fac= fac_inter[[19]] ,only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(FacCompSparseSp5$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(FacCompSparseSp5$success_comp,3)))
## Success rate for competition: 1
cat(sprintf("Success rate for facilitation: %s\n", round(FacCompSparseSp5$success_fac,3)))
## Success rate for facilitation: 0
data<-sim_data$FacCompSparseSp10
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjamcmp_facs10.rda",interact= (-1)*comp_inter[[20]] +fac_inter[[20]])
Rho_sign<-load_gjam(data,name="./gjam_models/gjamcmp_facs10.rda", interact= (-1)*comp_inter[[20]] +fac_inter[[20]])
g_metric_FacCompSparseSp10<-metrics_gjam(Rho_sign,comp=comp_inter[[20]],fac=fac_inter[[20]], only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_FacCompSparseSp10$success_env,3)))
## Success rate for non-interacting: 0.415
cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompSparseSp10$success_comp,3)))
## Success rate for competition: 1
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_FacCompSparseSp10$success_fac,3)))
## Success rate for facilitation: 0
data<-sim_data$FacCompSparseSp10
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hmcmp_facs10.rda" )
hm_conv(hm_mod)
Rho_sign<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact = (-1)*comp_inter[[20]] +fac_inter[[20]])
h_metric_FacCompSparseSp10<-metrics_hmsc(Rho_sign,comp=comp_inter[[20]],fac=fac_inter[[20]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompSparseSp10$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompSparseSp10$success_comp,3)))
## Success rate for competition: 0
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompSparseSp10$success_fac,3)))
## Success rate for facilitation: 0
data<-sim_data$FacCompSparseSp20
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjamcmp_facs20.rda",interact= (-1)*comp_inter[[21]]+fac_inter[[21]])
Rho_sign<-load_gjam(data,name="./gjam_models/gjamcmp_facs20.rda", interact= (-1)*comp_inter[[21]]+fac_inter[[21]])
g_metric_cmp_facs20<-metrics_gjam(Rho_sign,comp=comp_inter[[21]], fac=fac_inter[[21]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_cmp_facs20$success_env,3)))
## Success rate for non-interacting: 0.363
cat(sprintf("Success rate for competition: %s\n", round(g_metric_cmp_facs20$success_comp,3)))
## Success rate for competition: 0.5
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_cmp_facs20$success_fac,3)))
## Success rate for facilitation: 0.75
data<-sim_data$FacCompSparseSp20
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hmcmp_facs20.rda" )
hm_conv(hm_mod)
Rho_sign<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact = (-1)*comp_inter[[21]] +fac_inter[[21]])
h_metric_cmp_facs20<-metrics_hmsc(Rho_sign,comp=comp_inter[[21]],fac=fac_inter[[21]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_cmp_facs20$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(h_metric_cmp_facs20$success_comp,3)))
## Success rate for competition: 0
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_cmp_facs20$success_fac,3)))
## Success rate for facilitation: 0